#!/usr/bin/perl -w
use strict;
use Getopt::Long;
use lib "/net/cpp-group/Leo/bin";
use parse_fasta;
#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888
my $usage = <<"USAGE";

USAGE:

    fa_extract_domain_regions_using_hmm.pl
                        --hmm       an.HMM
                        --sequences seq.fa
                        --extra EXTRA_RESIDUES
                         [--regex "(\\S+)"]

                    > result.fa

        This script accepts sequences in the FASTA format and extracts
        the specified domain region and EXTRA_RESIDUES number of extra
		residues on each side.
USAGE

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory


# optional parameters
my $help = undef;
my $hmm = undef;
my $regex = "(\\S+)";
my $extra = 5;
my $seq_file = undef;
GetOptions('help' => \$help,
		   'hmm=s' => \$hmm,
		   'regex=s' => \$regex,
		   'extra=s' => \$extra,
		   'sequences=s' => \$seq_file);

die $usage if ($help);
die $usage unless ($hmm && $seq_file);

my %sequences;
sub save_seq_callback($$$$$$)
{
    my ($acc, $seq, $save_sequences) = @_[0, 1, 5];

	my $sequence = join "", @$seq;
	$sequence =~ tr/BJOUZ//d;
	if ($acc =~ /$regex/)
	{
		$sequences{$1} = $sequence;
	}
	else
	{
		$sequences{$acc} = $sequence;
	}
}
{
	open SEQ, $seq_file or die "Error:\n\tCould not open $seq_file\n$!";
	my $cnt_seq =parse_fasta::parse_sequences(  *SEQ,
												\&save_seq_callback,
												undef);
	close SEQ;
	print STDERR "$cnt_seq\tsequences parsed.\n";

}
open HMMOUTPUT, "/net/cpp-group/tools/hmmer/hmmsearch $hmm $seq_file |" or die;

my $cnt_hmm_seq = 0;
my $domain_name;
my $seq;
PARSE_DOMAIN: while (<HMMOUTPUT>)
{
	#next unless /Parsed for domains/;
    #dog_r1: domain 1 of 1, from 113 to 174: score 9.3, E = 0.00012
	if(/(.*):\s+domain\s+(\d+)\s+of\s+(\d+),\s+from\s+(\d+)\s+to\s+(\d+):\s+score\s+([-+]?\d+\.?\d?[eE]?[-+]?\d+),\s+E\s+=\s+([-+]?\d+\.?\d?[eE]?[-+]?\d+)/)
	{
		die unless (exists $sequences{$1});
		my $start = $4 - 1;
		my $len = $5 + 1 - $start;

		my $start_ext = $start - $extra;
		my $len_ext = $len + $extra + $extra;
		if ($start_ext < 0)
		{
			$len_ext += $start_ext;
			$start_ext = 0;
		}
		$cnt_hmm_seq++;
		$domain_name = $1;
		print ">$domain_name\n".substr($sequences{$domain_name}, $start_ext, $len_ext)."\n";
		#print STDERR "$domain_name\t\t$start_ext\t$len_ext\t$start\t$len\n";
		$seq = '';
		while (<HMMOUTPUT>)
		{
			if(/(.*):\s+domain\s+(\d+)\s+of\s+(\d+),\s+from\s+(\d+)\s+to\s+(\d+):\s+score\s+([-+]?\d+\.?\d?[eE]?[-+]?\d+),\s+E\s+=\s+([-+]?\d+\.?\d?[eE]?[-+]?\d+)/)
			{
				$seq =~ tr/a-z//d;
#				print ">$domain_name\n$seq\n";
				redo PARSE_DOMAIN;
			}
			if (/\s+$domain_name\s+\d+\s+([A-Za-z\-]+)\s+\d+/)
			{
				$seq .= $1;
			}
		}
		last;
	}
	#print STDERR $_;
	#exit;
}
if (defined $domain_name)
{
#	print ">$domain_name\n$seq\n";
}

print STDERR "$cnt_hmm_seq\tsequences with HMM hits.\n";
